% 物理和模型参数
rho = 1050; % 猪肝密度 (kg/m^3)
c = 3500; % 比热容 (J/kg·K)
kappa = 0.5; % 热导率 (W/m·K)
P = 30; % 功率 (W)
a = 0.01; % 长轴 (m)
b = 0.008; % 短轴 (m)
c_depth = 0.03; % 深度 (m)
V = (4/3) * pi * a * b * c_depth; % 椭球体积 (m^3)

% 网格参数
r_max = 0.05; % 最大r (m)
z_max = 0.05; % 最大z (m)
dr = 0.0001; % r的步长 (m)
dz = 0.0001; % z的步长 (m)
dt = 0.01; % 时间步长 (s)
t_max = 15; % 最大时间 (s)

% 网格定义
r = -r_max:dr:r_max;
z = -z_max:dz:z_max;
t = 0:dt:t_max;
Nr = length(r);
Nz = length(z);
Nt = length(t);

% 初始化温度场
T = 25 * ones(Nr, Nz, Nt); % 所有初始温度设为25℃

% 设置初始温度条件
P1 = [find(abs(r-0.002322) == min(abs(r-0.002322))), find(z==0)];
P2 = [find(r==0), find(abs(z-0.005288) == min(abs(z-0.005288)))];
P3 = [find(abs(r-0.003) == min(abs(r-0.003))), find(z==0)];
P4 = [find(abs(r-0.1) == min(abs(r-0.1))), find(z==0)];
P5 = [find(r==0), find(abs(z-0.1) == min(abs(z-0.1)))];

T(P1(1), P1(2), 1) = 55;
T(P2(1), P2(2), 1) = 55;
T(P3(1), P3(2), 1) = 100;
T(P4(1), P4(2), 1) = 25;
T(P5(1), P5(2), 1) = 25;

% 求解方程
for k = 1:Nt-1
    for i = 2:Nr-1
        for j = 2:Nz-1
            % Crank-Nicolson Method
            T(i, j, k+1) = T(i, j, k) + dt * ...
                (kappa/(rho*c)) * ((T(i+1, j, k) - 2*T(i, j, k) + T(i-1, j, k)) / dr^2 + ...
                (T(i, j+1, k) - 2*T(i, j, k) + T(i, j-1, k)) / dz^2);
        end
    end
end

time_step=700;
T_at_time1 = T(:, :, time_step);
T_adjusted = T_at_time1 - 25;

% 归一化到25到100的范围
T_min = 0; % 因为已经减去了25，所以最小值从0开始
T_max = 75; % 原最大值100，减去25后为75

% 归一化处理
T_normalized =((T_adjusted - T_min) / (T_max - T_min)) * 100;
T_at_time1 = T_normalized
% 现在T_normalized中的数据已经被归一化到25到100的范围内
figure;
r_corrected = r * (2.27/0.0075)* (2.7/2);  % 横坐标 (r) 进行校正
z_corrected = z * (2.27/0.0075)* (2.7/2);  % 纵坐标 (z) 进行校正

new_image = zeros(1001 + 228, 1001);

for i = 1:228
    new_image(i:i+size(T_at_time1, 1)-1, :) = new_image(i:i+size(T_at_time1, 1)-1, :) + T_at_time1;

end
%new_image = new_image / max(new_image(:)) * 100;
imagesc(r_corrected, z_corrected,  new_image);  % 绘制二维温度场
% caxis([30 40]);

axis xy;  % 设置坐标轴方向正确
xlabel('r (mm)');
ylabel('z (mm)');
title(['Temperature Distribution at t = ', num2str(time_step*0.01), ' s']);
colorbar;  % 显示颜色条

% 设置 z 轴的范围
%zlim([min(T_at_time(:)), max(T_at_time(:))]);  % 设置z轴的范围
% zlim([0, 100]);  % 归一化后最大温度为100
set(gca, 'ztick', []);  % 隐藏 z 轴刻度
set(gca, 'zticklabel', []);  % 隐藏 z 轴刻度标签

axis equal;  % 保证坐标轴比例相同
% T_at_time = T(:,:,time_step);
% T_at_time1 = T(:, :,time_step);  % 归一化温度场
% 
% % 校正r和z坐标
% r_corrected = r * (2.27/0.0075) * (2.7/2);
% z_corrected = z * (2.27/0.0075) * (2.7/2);
% 
% % 绘制温度场
% figure;
% imagesc(r_corrected, z_corrected, T_at_time1);  % 绘制二维温度场
% axis xy;  % 设置坐标轴方向正确
% xlabel('r (mm)');
% ylabel('z (mm)');
% title(['Temperature Distribution at t = ', num2str(time_step), ' s']);
% colorbar;  % 显示颜色条
% 
% % 移除z轴，隐藏z轴标签
% zlim([min(T_at_time(:)), max(T_at_time(:))]);  % 设置z轴的范围
% set(gca, 'ztick', []);  % 隐藏z轴刻度
% set(gca, 'zticklabel', []);  % 隐藏z轴刻度标签
% 
% axis equal;  % 保证坐标轴比例相同
% % 可视化最终温度分布
% [X, Z] = meshgrid(r, z);
% surf(X, Z, squeeze(T(:,:,end))');
% xlabel('r (m)');
% ylabel('z (m)');
% zlabel('Temperature (°C)');
% title('Final Temperature Distribution');
